Dark solitons near the Mott-insulator— superfluid phase transition 
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Dark solitons of ultracold bosons in the vicinity of the Mott-insulator-superfiuid phase transition 
are studied. Making use of the Gutzwiller ansatz we have found antisymmetric eigenstates corre- 
sponding to standing solitons, as well as propagating solitons created by phase imprinting. The 
properties of both types of solitons are studies and in particular it is demonstrate that quantum 
fluctuations greatly affect the characteristics of the solitons. 
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Introduction. - Ultracold atomic gases provide a perfect 
playground to study nonlinear atom optics, and nonlin- 
ear structures and textures, such as solitons [H,B. These 
studies have led to the observations of dark [H, H[ and 
bright H} solitons in trapped Bose- Einstein condensates, 
or more recently, bright solitons stabilized by the pres- 
ence of dark ones @ . The analogy to nonlinear optics 
has triggered theoretical interest in discrete (lattice) soli- 
tons [3, II , and has led to the seminal observations of gap 
solitons, i.e., lattice solitons with repulsive interactions, 
but with an appropriate dispersion management Q. 

While most of the studies of solitons were concen- 
trated on their classical aspects, more recently, consid- 
erable interest has been devoted to the role of thermal 
noise [1, [l(| , quantum properties of solitons and the role 
of quantum fluctuations. The latter may cause filling up 
of the dark soliton core in the quantum detection process, 
as was shown using the Bogoliubov-de Gennes (BdG) 
equations 11]. The BdG method was also employed to 
study the stability of solitons [lH, excitations caused by 
the trap opening [l3j . and entanglement generation in 
collisions of two bright solitons [14| . A noisy version of 
the standing bright solitons was studied using the exact 
diagonalization and quantum Monte Carlo method [ToT ]. 
Bright soliton in ID were considered in Ref. [l3j], where 
exact Lieb-Liniger solutions were used to calculate the 
internal correlation function of the particles positions. 
Making use of the discrete nonlinear Schrodinger equa- 
tion (DNSE) , and the time-evolving block decimation al- 
gorithm 1 1 61 ] it was demonstrated that quantum effects 
lead the soliton to fill in, and that soliton collisions be- 
come inelastic (l7j ]. 

All the previous studies of lattice solitons were done 
deep in the superfluid phase. However, near the phase 
boundary between superfluid (SF) and Mott insulator 
(MI) [lq . the propagation of matter waves becomes 
strongly suppressed due to the enhancement of quan- 
tum fluctuations. In general, the properties of solitons 
in this regime are expected to be very different from the 
ones far in the SF regime. The aim of the present work 
is to study lattice solitons near the MI-SF phase tran- 



sition. This is done employing the position-dependent 
Gutzwiller ansatz which gives a satisfactory description 
of the quantum phases of inhomogeneous bosonic sys- 
tems (13 . [20l | as well as their dynamical behavior [2l| . 
The method was recently utilized to study vortices in 
the vicinity of the MI-SF phase transition [221 ]. 

The model.- We consider a system of ultracold inter- 
acting bosons in a d-dimensional lattice described by the 
Bose-Hubbard Hamiltonian 
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where J is the tunneling matrix element, U is the on- 
site atom-atom interaction energy, and \x the chemical 
potential. Throughout the paper, we will be dealing with 
repulsive interaction, i.e., U > 0. The annihilation and 
creation operators at site i, cti and a!, obey the bosonic 
commutation relations. 

Our analysis employs the Gutzwiller ansatz. Thereby, 
eigenstates of the Hamiltonian ([TJ) are taken as products 

of local states |$) = riil s i)> l s i) = X^=0 c m| n )i- Here > 
\n)i is the Fock state with n atoms at site i. We are in- 
terested in the solutions which has a position-dependence 
only in one spatial dimension. Then the equations of mo- 
tion [H, H3, (H take the form 
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yfn is the condensate order 



where ipi = 
parameter. 

Standing solitons.- Stationary solutions of Eqs. @ 
are worked out for a lattice with the finite number of 
sites L and with the boundary conditions co, (l = ci,„, 
CL+i.n = cl.u, for any n. In the ground state, the co- 
efficients Cm and the order parameter tpi do not depend 
on the site index i. The boundaries between the MI and 
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FIG. 1: Mean number of atoms (hi) (a) and (c), and mean- 
field order parameter tpi (b) and (d). The scaled chemical 
potential n/U = 1.2 and the tunneling rates 2dJ/U: 0.7 (i), 
0.5 (ii), 0.3 (iii), 0.15 (iv), and 0.05 (v). 
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FIG. 2: Dark areas bounded by the solid lines show the 
first three MI zones. Gray areas between solid and dotted 
lines indicate the regions where the off-site solitons have a 
local maximum of the mean occupation number at the middle 
lattice sites. The light gray areas between solid and dashed 
lines depict the same for on-site solitons. In the rest part 
of the diagram the mean particle number takes the minimal 
value at the middle site(s). 



SF are determined by 2dJ c /U = (n — ijl/U)(h/U — n + 
1)/(1 + n/U), where no is the smallest integer greater 
than n/U. For J < J c , the ground state is MI, where 
Cm = S n , no exp(-iw n t), huj n = Un(n-l)/2- (in, tpi = 0. 

In the present work, we are interested in the low-energy 
excited eigenstates, where the order parameters tpi are 
anti-symmetric with respect to the middle point of the 
lattice. These are the kink states which can be treated 
as standing dark solitons. In the Gutzwiller ansatz, the 
excited states of the MI are products of local Fock states, 
where the occupation numbers rii can be different from 
no. As a consequence, all ipi vanish again and the soli- 
ton solutions are lacking for the values of the parameters 
within this regime. 

In the SF, the situation is different and one has to 
distinguish between two cases: when the middle point is 
on the lattice site (onsite modes) and in the middle of two 
neighboring sites (off-site modes). The two modes have 
different energies, and the difference defines the Peierls- 
Nabarro barrier (24[. In the present work, we will be 
dealing only with situations where the presence of the 
barrier is not relevant. 

Typical examples of standing onsite as welll as off-site 
soliton modes are displayed in Fig. [TJ The mean oc- 
cupation numbers {hi) are shown in (a) and (c), while 
(b) and (d) give the associated mean-field order pa- 
rameters ipi. The individual curves correspond to dif- 
ferent tunneling rates J. For the considered chemical 
potential, /i/U = 1.2, the MI-SF transition occurs at 
2dJ/U » 0.0727. Above this value, in the SF regime, 
solitons are indeed found as depicted in Fig.QJb) and (d). 
In the numerical calculations, n in Eq. @ is restricted by 
some finite N (c^jv+i = 0). The computations are per- 
formed using the imaginary time propagation technique. 



The size of the lattice L as well as the cut-off number of 
atoms N at each site where chosen large enough such that 
their influences on the eigenstates are negligible. The di- 
mension is d = 3 for these examples and throughout the 
paper. 

In the present work, we consider the parameters range 
near the MI-SF transition where quantum fluctuations 
are especially important. An outcome of this is that 
the mean occupation numbers (rii) are not equal to the 
mean number of condensed atoms \ipi\ as it would in the 
mean-field description valid deep in the SF phase. The 
relation between ipi and the expansion coefficients Cj >n is 
given after Eq. ([2]). For example, the non-zero occupa- 
tion number, despite tpi = 0, is clearly seen from Fig. [T] 
(a) and (b). For the onsite mode, when ipi vanishes at the 
middle site, the state of the corresponding site is a Fock 
state. In the example of Fig. [IJ this is the Fock state with 
n = 2. On the other hand, for the off-site mode we have 
ipi 7^ for all i's in the SF regime and thereby the onsite 
states are never pure Fock states. Far in the SF regime, 
the occupation numbers attain a global minimum at the 
middle lattice site(s). Surprisingly, by approaching the 
MI boundaries the global minimum turns into two lo- 
cal minima at the neighboring sites. Even closer to the 
critical point, a global maximum is formed at the mid- 
dle lattice site(s). We distinguish these localized modes 
from the ones with a single global minimum at the mid- 
dle site(s). Fig. [2] depicts a (/x, J) diagram identifying 
the various types of solutions. The anomalous regions 
where the mean occupation numbers have local maxima 
are located near the lower branches of the MI- lobes corre- 
sponding to the hole excitations [25| which support dark 
solitons. 

Propagating solitons created by phase imprinting. - Ex- 
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FIG. 3: Time evolution of the mean occupation numbers {hi) 
(a) and the order parameters \^i\ 2 (b) after phase imprinting 
with A<f) — ir, kmp = 2. The parameters are fi/U = 1.2, 
2dJ/U = 0.3, giving similar evolution for both {hi) and lipif- 
The spatial region (vertical axis) contains 80 lattice sites, and 
the dimensionless time t — tU /h. 



perimentally, dark (or gray) solitons are typically created 
via a phase-imprinting method [3j. Initially (t — 0) the 
system of atoms is assumed to be in its ground state. 
During a short time U mp one applies a spatially depen- 
dent potential on top of the lattice. In the Bose-Hubbard 
Hamiltonian, it is described by the term ejdjcij. If the 
time Ump is much shorter than other characteristic time 
scales, from Eqs. © we get that it induces a shift in the 
phase of the atomic states 



QniUmp) = Cfa(0)exp ( „i £ ^""P , 



ipi(Ump) = ipi(0)exp l-i 
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For the creation of dark solitons we choose a hyperbolic 
tangent imprinting potential, such that 



A0 
~2~ 



1 + tanh 



I 



0.45Z 



(7)1 p 



(4) 



Here, Ump is the width of the interval around I = 0, 
where grows from 0.1 to 0.9 [26(. Apart from the 

moving gray soliton, the phase imprinting also induces a 
density wave propagating in the opposite direction to the 
soliton HIHl. 

Figs. [31 HI [5] present the time evolution of the occu- 
pation numbers (hi) and squared absolute values of the 
mean- field order parameters IV'il 2 for two different values 
of J and fi. Scaling fi and J by U, we work with a dimen- 
sionless time t — tU/ft. In Fig. [3J J is taken relatively 
large but still in the regime where the on-site modes have 
local maxima of {hi) (light gray region in Fig. [21). In this 
case, the outcomes of quantum fluctuations are small and 
the dips and the maxima of {hi) coincide with those of 



FIG. 4: Same as in Fig. H but with 2dJ/U = 0.15. In this 
regime, the directions of propagating maximum and minimum 
of \ipi\ 2 and (rij) seen in Fig. [3] are reversed, r is the scaled 
dimensionless time. 



\ipi\ 2 - We note that the dips propagate slower than the 
maxima. This type of dynamics is always recovered in 
simulations based on the DNSE. 

In Fig. HI the value of /i is the same as in Fig. [3] but 
J is smaller such that we enter into the gray region of 
Fig. [21 The overall structure of {hi) remains the same as 
in the example shown in Fig. [3j but here \ipi\ 2 behaves 
differently. More precisely, \ipi\ 2 shows a local maximum 
coinciding with the propagating dip of the occupation 
numbers, and a local minimum where instead {hi) has 
a local maximum. Furthermore, the dip of (rij) propa- 
gates faster than its maximum moving in the opposite 
direction. This anomalous behavior, which is found only 
within the gray regions of Fig. [21 cannot be described by 
the DNSE. Closer to the boundary of the gray region as 
in the example depicted in Fig. [5l |"0i| 2 can become os- 
cillating and spreading around the imprinting site, while 
{hi) still shows similar dynamics as in Figs. [3J and [U In 
the remaining white regions of Fig. [21 the time evolution 
is qualitatively the same as in Fig. [3j 

We also performed simulations with other values of 
limp and did not find any strong influence on the propaga- 
tion velocity. Larger values of li mp result in broader prop- 
agating modes and their shapes become more smooth. 
The interference pattern between the modes propagat- 
ing in the opposite directions, visible in the center of 
Figs. El [H O becomes suppressed for larger k mp . 

Finally, the soliton velocity v so i as well as the velocity 
of the density wave Vdens as functions of A(j) are shown in 
Fig. [6] They are calculated making use of a linear fit of 
the corresponding minimum and maximum of {rii) for the 
dimensionless time r = tU/fi > 10. Although the depen- 
dences are very weak, it is found that the velocities are 
not monotonic showing a maximum for A</> ps 7r/4. This 
behaviour is different from what was found in the simu- 
lations based on the Gross-Pitaevskii equation in contin- 
uum in a harmonic trap Q or in a periodic potential [l^ | . 
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FIG. 5: Same as in Figs. |3] and 0] but with n/U = 1.4, 
2dJ/U = 0.15. For these parameters, \ipi\ 2 oscillates and 
spreads near the imprinting site. Time r is dimensionless. 
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FIG. 6: Dimensionless velocity of the soliton (solid line) 
and density (dashed line) modes as a function of A(j>. Here 
2dJ/U = 0.3, fj,/U = 1.2, and l imp = 2. 



With the decrease of A(j> the propagation velocity slightly 
increases (Fig. |6]). On the other hand, the depth of the 
propagating dip as well as the height of the maxima de- 
creases. The same was found in continuum model [3[. 

Conclusion.- We have investigated dark solitons of 
bosons in optical lattices at zero temperature. Using 
the Gutzwiller ansatz, we found kink stationary states 
where the condensate order parameter ipi is antisymmet- 
ric function with respect to some spatial point. However, 
in certain regions close to the MI-SF transition which 
are shown in Fig. [21 the corresponding ni does not have 
a global minimum at the point where ipi vanishes. This 
anomalous behavior shows up near the phase boundary 
where quantum fluctuations play an important role. 

The real-time dynamics of the propagating dark soli- 
tons created by the phase imprinting is studied as well. 
In the MI phase, the solitons cannot be created. This 
can be done only in the SF phase, where there are al- 
ways global minima and maxima of n.; propagating in 
the opposite directions and these directions are always 
the same. The behavior of the condensate can be dif- 
ferent and this happens in the anomalous regions of the 
off-site stationary modes shown as gray areas in Fig. [2] 
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